Phase transitions in spin-orbital models with spin-space anisotropics for 
iron-pnictides: A study through Monte Carlo simulations 
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The common phase diagrams of superconducting iron pnictides show interesting material specifici- 
ties in the structural and magnetic phase transitions. In some cases the two transitions are separate 
and second order, while in others they appear to happen concomitantly as a single first order transi- 
tion. We explore these differences using Monte Carlo simulations of a two-dimensional Hamiltonian 
with coupled Heisenberg-spin and Ising-orbital degrees of freedom. In this spin-orbital model, the 
finite-temperature orbital-ordering transition results in a tetragonal-to-orthorhombic symmetry re- 
duction and is associated with the structural transition in the iron-pnictide materials. With a zero 
or very small spin space anisotropy, the magnetic transition separates from the orbital one in tem- 
perature, and the orbital transition is found to be in the Ising universality class. With increasing 
anisotropy, the two transitions rapidly merge together and tend to become weakly first order. We 
also study the case of a single-ion anisotropy and propose that the preferred spin-orientation along 
the antiferromagnetic direction in these materials is driven by orbital order. 

PACS numbers: 74.70.Xa,75.25.Dk,75.40.Cx,75.40.Mg 



INTRODUCTION 

Several parent phases of the superconducting iron- 
pnictide materials show an interesting interplay of 
structural, magnetic and orbital degrees of free- 
dom. [l|43| These quasi two-dimensional (2D) mate- 
rials share a similar phase diagram, where a tetrago- 
nal paramagnetic phase at high temperatures transitions 
into an orthorhombic, antiferromagnetic phase at low 
temperatures. [3-0] The square-lattice of iron atoms de- 
velop magnetic order at wavevector (7r,0), which corre- 
sponds to antiferromagnetic (AFM) alignment of spins 
along one of the nearest-neighbor direction (x) and ferro- 
magnetic (FM) alignment along the other (y) . [lClMl9l | The 
orientation of the ordered spin moments is tied to antifer- 
romagnetism and points along the AFM direction. 13- 19j| 

While lattice distortions are typically quite small in 
iron pnictides, the observed spin- wave spectra from neu- 
tron scattering suggest a robust, possibly sign changing 
anisotropy in the exchange constants along the x and y 
directions. [20l - |22| Various transport, optical, and spec- 
troscopic measurements also show substantial emergent 
anisotropics in the 2D xy plane. [23r,-32j In particular, 
an orbital polarization associated with the occupation 
of dxz and dyz orbitals has been observed. 33-36| These 
anisotropics in some cases can persist up to high temper- 
atures and have been identified with long and sometimes 
short range Ising- nematic order. 37, 11] 

Despite the above similarities, there are also substan- 
tial material-specific differences. The parent compounds 
of the 1111 family (RFeAsO, with R a rare earth element) 



of iron pnictides undergo two separate second order phase 
transitions, where the structural transition is followed by 
a magnetic transition at a lower temperature. Q On the 
other hand, in the 122 family (AFe2As2, with A an alka- 
line earth element) the two transitions appear to occur 
at the same temperature. 

More recent measurements revealed that in the un- 
doped BaFe2As2, the structural and magnetic transitions 
are slightly separated by less than one Kelvin. [Mii In 
that case, the structural transition starts as second order, 
and at a slightly lower temperature there is a first order 
jump in the lattice distortion with a concomitant first 
order magnetic transition. This feature is not generic to 
all the 122 family of iron pnictides. In particular, there is 
strong evidence showing a largely first order phase tran- 
sition in CaFe2As2 and SrFe2As2, where the structural 
and magnetic phase transitions coincide. 42h44 

There have been many proposals for the mechanism 
driving these transitions. These include, (i) emer gent 
Ising nematic orders in frustrated spin systems, [37|, |38 . 

-[SSj (u) orbital order J5j-l63j (iii) coupling to lattice 
degrees of freedom, 6J, |65| and (iv) symmetry break- 



ing associated with fermi-surface effects in an itinerant 
system. 65 -7^ On symmetry grounds one cannot distin- 
guish between different pictures, since the different de- 
grees of freedom lead to same broken symmetries and 
they are all present to some extent and coupled to each 
other. Thus, detailed quantitative studies are important 
to establish the role played by different mechanisms. 

In this paper, we wish to study the scenario where 
orbital order is the primary driving mechanism for the 
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finite temperature transitions. We investigate the prop- 
erties of a spin-orbital model, where the spin and orbital 
degrees of freedom are coupled by a Kugel-Khomskii like 
mechanism. 71| In the model, the local orbital occupa- 
tion modulates the spin exchange constants. Once the 
orbitals are ordered, coUinear antiferromagnetism can de- 
velop and anisotropic exchange constants in the x and y 
directions Jix ^ Jiy result. However, we note that a 
model containing only effective Heisenberg spin interac- 
tions JijSi ■ Sj) is still isotropic in spin space, since 
the energy of the system does not depend on the direction 
of magnetization with respect to the crystal axes. 

A simple mean-field treatment of our spin-orbital 
model suggests that the spin and orbital orderings oc- 
cur simultaneously as a single phase transition, which 
can be first or second order depending on the ex- 
change couplings. However, such a treatment neglects 
long-wavelength fluctuations which can drive the spin- 
ordering temperature to zero. To study the effects of 
fluctuations, we employ large scale Monte Carlo simu- 
lations by treating the spin and orbital variables clas- 
sically, which should be sufficient for finite-temperature 
phase transitions. The Monte Carlo results indicate that 
if spin rotational invariance is preserved, at finite temper- 
atures there is only one orbital ordering transition which 
belongs to the 2D Ising universality class. In this case, 
long-range spin order only occurs at T = 0, in accord 
with the Mermin- Wagner theorem. However, a small spin 
space anisotropy ( 5%) will bring the magnetic transi- 
tion temperature up to the orbital one. With increasing 
anisotropy the coupled spin-orbital transition tends to 
become first order. These results are reminiscent of the 
observed behaviors of different families of iron pnictides. 

Our study neglects three-dimensional (3D) couplings, 
studied for example in Refs. I38l.l5ll andlSSl 3D couplings 
have a similar effect as spin space anisotropics in that 
they both can result in a finite-temperature magnetic 
transition. However, in general they will lead to different 
universality classes for the transitions. While 3D cou- 
plings could be more important in some materials (for 
example within the 122 family), spin space anisotropy 
may be more important in others. In some materials the 
magnetization has been reported to obey 2D Ising uni- 
versality behavior. I76i| Even if the ultimate transition is 
weakly first order in these materials, the reported fluc- 
tuations appear more 2D. This provides a motivation for 
our choice of anisotropy over 3D couplings. 

Studying spin space anisotropy also allows us to ad- 
dress the orientation of the ordered moments. In a quasi- 
two dimensional material, one would expect the uniax- 
ial anisotropy to point out of the plane and the spins 
should be equally likely to point along any direction in 
the plane. However, this can change with orbital polar- 
ization. In transition metal compounds, ligand crystal- 
field splitting can lift the degeneracy of the transition 
metal id orbitals. In this case, the orbital moments are 



usually quenched and there may be no preferred spin 
directions. However, relativistic spin-orbit coupling can 
induce a non-zero orbital angular momentum, which ac- 
companied by an orbital polarization (such as a prefer- 
ential occupation of dxz over dyz orbitals) can lead to 
a single-ion anisotropy term and an anisotropic ij-factor 
in the xy plane. In this case, excess population of d^z 
orbitals can favor spins pointing along the x-axis, while 
excess population of dyz orbitals will favor spins pointing 
along the y direction. We propose that in iron pnictides 
the single-ion anisotropy term in the xy plane is related 
to orbital order and since it is also tied to AFM it leads 
to spins pointing along the AFM direction. 

The outline of the paper is as follows. In section II 
we introduce our model and develop a simple mean field 
theory. In section HI the Monte Carlo results are dis- 
cussed. In section IV we discuss the implication of our 
study for the iron pnictide materials and in section V we 
summarize our work. Details of the Monte Carlo method 
are presented in the appendix. 



II. SPIN-ORBITAL MODEL 

The spin-orbital model is given by the Hamiltonian, 

H =^^{Jin^ni+x - JF)Si ■ Si+x (1) 

■i 

+ ^(Ji(l - ni){l - n,,+y) - Jf)Si ■ Si+y 

i 

+ J2Si ■ Sj . 

Here Si are classical Heisenberg spins on a square-lattice. 
Hi are classical Ising variables that take values or 1, 
and << i,j >> signifies summing on next nearest neigh- 
bor pairs of the square-lattice. In a classical system the 
spin magnitude S is not important in determining the 
phase transitions, and for our discussion we set S* = 1. 
Physically, the variables Ui represent the preferential oc- 
cupation of dxz [ni — 1) or dyz {ui — 0) orbitals. The 
model has tetragonal symmetry. However, Ui — 1 (occu- 
pation of dxz orbitals) favors AFM order along the x-axis, 
whereas = (occupation of dyz orbitals) favors AFM 
order along the y-axis. We have added an orbital inde- 
pendent FM nearest neighbor interaction Jp and used 
J2 ~ 0.4Ji and Jp = 1/6 Ji. This set of parameters cor- 
responds to the neutron scattering observation that the 
spin-wave spectra is better fit with an AFM exchange 
along one direction and a weak FM exchange along the 
other. [tI] The latter could arise from double exchange 57 
or from the orbital geometries. |55i] But, its sign or mag- 
nitude is not crucial for the phase transitions we report 
here. Spin-space anisotropics will be introduced later 
when we discuss the Monte Carlo simulations. 
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The ground state of this model breaks tetragonal sym- 
metry. It has a ferro-orbital order, all = 1 or = 0, 
corresponding to nearest neighbor exchanges which are 
AFM along one axis and FM along the other. The ground 
state has (7r,0) spin order when rii = 1 and (0,7r) spin 
order when = 0. 

We note here that in iron pnictides, the low tempera- 
ture orbital polarization is found to be incomplete, where 



the occupation number is not strictly one or zero. [33143 



A partial orbital polarization can result from the itiner- 
ant electron degrees of freedom, or from quantum fluctu- 
ations in the orbital variables. The role of orbital order 
in driving the structural and magnetic transitions of iron 
pnictides indeed has been discussed based on an itin- 
erant electron perspective using multi-orbital Hubbard 
Hamiltonians. [Hi-lil Our approach of studying a Kugel- 
Khomskii type spin-orbital model can be viewed as the 
strong coupling limit of such Hamiltonians. While we 
leave out the charge degrees of freedom which are impor- 
tant in describing for example transport properties, our 
model should still capture the key physics of magnetism 
and finite temperature phase transitions. 

Below we first develop a mean-field theory for the 
phase transitions of the spin-orbital model under consid- 
eration. We set rii = {l + ai)/2 and assume a mean-field 
Hamiltonian of the form 



n 



MF 



^-Y,B'^S,i-Y,B2Sa-hY,a, (2) 



where the first sum runs over sublattice one, the second 
over sublattice two, and the third over all the spins in 
the lattice. Bi and B2 are the staggered fields on the 
two sublattices and /i is a field that couples to orbital 
order. Focusing on the (tt, 0) order, we let m = {Si) and 
n = (cr,) > 0. We find for i = 1, 2 

B,^2m{Jin + 2J2), /i = Jlm^ (3) 

leading to the mean-field equations, 

m = L(2l3m{Jin + 2J2)), (4) 

and 



tanh(/3JiTO^). 



(5) 



Here, L{x) is the Langevin function cothx — 1/x. These 
equations lead to a simultaneous transition and an orbital 
ordered AFM phase. It is a second order phase transition, 
with a transition temperature of 4J2/3, provided J2 > 
constant Ji. The transition becomes first order when Ji 
exceeds J2 (the case of interest in the pnictides) . 

While mean-field theory can not be quantitatively valid 
because of the divergent infrared fluctuations in the spin 
variable, which push the spin ordering transition tem- 
perature to zero, we will see that the mean-field results 
correctly capture the following physics: 



1. Non-zero magnetic order produces an ordering field 
for the orbital degrees of freedom. Hence, whenever there 
is magnetic order present, orbitals symmetry will also be 
broken. Thus, orbital transition can not happen below 
the magnetic ordering transition. 

2. Without some order of the magnetic degrees of free- 
dom, the orbitals do not interact. Actually, orbital cou- 
plings depend on short-range magnetic order not long- 
range magnetic order. This is not allowed for in the 
mean-field theory but will become clear from our later 
discussion of the Monte Carlo simulations. Thus, the 
two transitions are always going to be close in temper- 
ature, unless the magnetic transition is pushed signifi- 
cantly below the mean-field transition temperature due 
to additional fluctuations. 

3. The orbital ordering temperature is not significantly 
depressed by the fluctuations of the spins and our mean 
field theory provides a fairly good prediction of the tran- 
sition temperature. 

4. We will see in the Monte Carlo simulations that 
the main role of the long-wavelength spin fluctuations 
is to decouple the spin and orbital transitions. The spin 
transition temperature is pushed to zero in the absence of 
spin space anisotropy. In this case, the orbital transition 
becomes Ising like and second order. 

5. With significant anisotropy, both the spin and or- 
bital transition temperatures rapidly approach the mean- 
field values, and the transition has a tendency to become 
first order for Ji > J2. 



III. RESULTS OF MONTE CARLO 
SIMULATIONS 

In this section, we present the results of Monte Carlo 
simulations with and without spin-space anisotropics. 
The details of the Monte Carlo methods as well as the 
quantities measured and the expected scaling behavior 
are discussed in the appendix. 



Isotropic Heisenberg Spins 

First, we consider the case of isotropic Heisenberg 
spins. The squares of spin and orbital order parame- 
ters obtained from the simulation are shown in Fig. [1] 
We know on general grounds that in a 2D system spin 
rotational symmetry can not be spontaneously broken 
at any finite temperature. However, this is not evident 
from the plot. The exponential growth of the spin-spin 
correlation length rapidly exceeds the size of the system 
and this creates the impression of long-range order at a 
finite temperature. One needs to carefully study the size 
dependence. The Binder ratios, defined in the appendix, 
prove useful for this purpose. 
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FIG. 1: Squares of spin and orbital order parameters as a 
function of temperature for the isotropic spin-orbital model 
on a 20 X 20 lattice. The vertical line shows the transition 
temperature Tc (measured in units of Ji), where the orbitals 
develop long-range order. 
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FIG. 3: Orbital Binder ratio as a function of temperature 
(measured in units of Ji) for different L x L lattices. In con- 
trast to the spin variables, the discrete orbital variables un- 
dergo a phase transition, developing long range order at finite 
temperatures. 
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FIG. 2: Spin Binder ratio as a function of temperature (mea- 
sured in units of Ji) for different L x L lattices. For the 
isotropic spin-orbital model there are no crossings in gs at 
any temperatures of our simulation. This is consistent with 
the theory that in a 2D system there is no long range spin 
order at finite temperatures. 



Figure m gives the spin binder ratios, gs, which show 
no crossings with system size down to the lowest mea- 
sured temperature, signifying absence of long range order 
at finite temperatures, in agreement with the Mermin- 
Wagner theorem. In contrast the orbital binder ratios 
gn, shown in Fig. [3l have clear crossings at finite tem- 
peratures and we can extract Tc by comparing different 
system sizes. We obtain T^/Jl 0.450 ± 0.001 for the 
isotropic spin-orbital model. 



In Fig. m we show the scaling plot for the orbital sus- 
ceptibility. The data collapse leads to estimates of criti- 
cal exponents ly = 1.01 ± .01 and 7 = 1.75 ± .02. These 
exponents are consistent with the 2D Ising universality 
class. Figure [S] shows a plot of the specific heat, which 
grows rapidly near Tc- It is consistent with a logarith- 
mic divergence but with an amplitude significantly larger 
than that in the pure 2D Ising model. The amplitude of 
the specific heat is not universal but is comparable for the 
Ising model on different 2D lattices. 75] It is considerably 
larger in our model, presumably as the Ising nematic vari- 
ables associated with the spins couple directly to the or- 
bitals and enhance the amplitude. We last note that the 
sharp peak in our specific heat clearly indicates a phase 
transition and its transition temperature T^,. Therefore, 
the Monte Carlo simulations of of our spin-orbital are 
less affected by finite-size effect compared to that in the 
frustrated square lattice J1-J2 model. 
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Exchange and single-ion anisotropy 

We next consider models with spin space anisotropy 
by generalizing the scalar product 



Dj — Oj Dj 



A[5f 5J + SfS^] 



(6) 



We study the system for several values of the Ising 
anisotropy parameter A. As long as A < 1, there is only 
Ising symmetry for the spins, and both spins and or- 
bitals can order at finite temperatures. The effect of the 
anisotropy on the orbital order is small, and the transi- 
tion temperature is raised gradually as A is reduced. In 
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FIG. 4; The scaling of the universal x versus reduced tem- 
perature \t\ = |(r — Tc)/Tc\ shows that critical exponents are 
consistent with the 2D Ising Model. 



FIG. 6: Spin Binder ratio as a function of temperature (mea- 
sured in units of Ji) for different LxL lattices. For anisotropic 
spins with either Ising or single ion anisotropy, finite spin or- 
dering is observed besides orbital order. 
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FIG. 5: Specific heat for the isotropic spin-orbital model un- 
der study. The sharp peak is consistent with a logarithmic 
divergence at Tc. 



contrast, one can see a dramatic difference in the Binder 
ratios for the spin variables. Comparing to isotropic 
spins, viz. A = 1, in Fig. [2l there are clear crossings 
in Fig. [S] We can extract for both order parameters 
using the Binder ratio. 

We summarize the extraction of Tc over a range of A in 
Fig. [71 When A is near (but not equal to) 1, we are in a 
regime where the spin transition temperature is non-zero 
but still separated from the orbital transition. However, 
when the anisotropy is small (in our case, |1 — A| < 0.1), 
computationally it is difRcult to distinguish the two tran- 
sitions. This shows that the spin transition tempera- 
ture grows very rapidly with increasing anisotropy and it 



FIG. 7: From Binder cumulant ratios systematic crossings are 
located for systems of size L and 2L. These are plotted versus 
inverse system length and extrapolated to get the thermody- 
namic Tc (measured in units of Ji). In all cases, the orbital 
crossings (solid) approach T from above and the spin cross- 
ings (dashed) approach it from below. For A = 1, spins order 
only at zero temperature. 

rapidly merges with the orbital transition. We note that 
with anisotropy the transition temperatures are within 
10% of the mean-field value of 0.53Ji. 

We now introduce a single ion anisotropy, which is tied 
to orbital order. 

N 

H^on = -DY,{n,Sf + (1 - n,)Sf) . (7) 

i 

In transition metal compounds, ligand crystal-field split- 
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FIG. 8: In the orbital configuration {rn — 1}, the spin-orbital 
model with single ion anisotropy has an AFM exchange along 
the X direction and favors the spin order collinear with this 
exchange. This is shown by plotting the 5*^ component from 
a typical spin configuration when in the {ui = 1} phase. The 
false color plot represents magnitude of the spin component 
along the x direction. 
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FIG. 9: For an isotropic ( A = 1 ) 24 x 24 system, we plot 
the orbital and spin binder ratios, g„ and gs- g-n behaves 
smoothly between its limiting values while gs develops what 
seems like a divergence at the transition temperature. The 
development of such an incipient divergence is an indicator of 
a first order transition. 



ting lifts the degeneracy of the transition metal 3d or- 
bitals, and the orbital angular moments are usually 
quenched. In this case, treating the relativistic spin-orbit 
coupling as a perturbation to the second order will result 
in a single-ion anisotropy term that reflects the underly- 
ing symmetry of the crystal. Therefore, an orthorhombic 
structural distortion or a net orbital polarization can lead 
to a single-ion anisotropy term closely tied to orbital or- 
der. This single ion anisotropy favors spin orientations 
along the AFM direction. One ordered configuration ob- 
served in the simulation is shown in Fig. [5] 

Similar to the case of Ising- anisotropy (Eq. ©), we 
have simulated these systems with different D values. 
For \D\ > 0.1, once again we see no separation of the two 
ordering transitions for spins and orbitals. This shows 
that, like the Ising anisotropy, the uniaxial anisotropy 
causes a rapid increase in the spin ordering temperature 
and it soon merges with the orbital order. 

In Fig. ini we show the binder ratios for orbitals and 
spins. Binder ratios for orbitals remain well behaved 
regardless of the anisotropy introduced in the models. 
However, there is a clear incipient divergence in gs, which 
is indicative of a first order transition. In general, we find 
that the spins have a greater tendency for a first order 
transition than the orbitals. The implications of these re- 
sults for the pnictides are discussed in the next section. 

V. DISCUSSION AND RELEVANCE TO THE 
IRON PNICTIDES 

In this section, we use the results of Monte Carlo 
simulations, mean field theory and general arguments 



about quasi-2D spin systems to develop an overall phase- 
diagram for coupled spin-orbital systems. We will then 
explore the applicability of the phase diagram to the iron- 
pnictide materials. The key issues of interest to us are 
whether there is a single transition or two separate tran- 
sitions, and whether each of the transitions is first or 
second order. 

The phase diagram of the spin-orbital model 

The mean-field theory gives a simultaneous spin and 
orbital transition, which could be first or second order de- 
pending on the exchange couplings. Monte Carlo simula- 
tions show that the spin and orbital transitions are prac- 
tically simultaneous unless the spin space anisotropy is 
very small. In the latter case, divergent long-wavelength 
fluctuations push the spin transition temperature to zero, 
whereas the orbital transition is not significantly affected 
by these fluctuations. The transition temperatures ob- 
served in the simulations are within a few percent of the 
mean-fleld value of 0.53 Ji when the anisotropy is large. 
As the anisotropy goes to zero, the orbital transition is 
reduced by less than 20%, whereas the spin transition is 
reduced all the way to zero. Even a 5% anisotropy causes 
a near simultaneous transition. 

On general grounds, one knows that in place of long- 
range order the correlation length in a 2D Heisenberg spin 
system stays finite but grows exponentially as expC/T 
as the temperature is lowered. This implies that if be- 
low some energy scale eg these divergent fiuctuations are 
cut off (due to for example spin space anisotropy or 3D 
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FIG. 10; Phenomenological phase diagram for the spin-orbital 
model. The exchange energy scale in the problem sets the 
transition temperature To for orbital order, which in turn 
drives the structural transition, eo is the energy scale below 
which long wavelength fluctuations are suppressed. There are 
two separate continuous orbital and magnetic transitions for 
e <SC £0 (shown as dotted and dashd lines respectively) and 
one simultaneous first order transition for e ^ eo (shown as a 
solid hue). Near the region e ~ eo (segment AB in the figure), 
the two transition temperatures can be very close and can be 
continuous or first order. In iron pnictides, e would refer to 
the larger of the spin anisotropics or 3D couphngs. 



coupling), it will lead to long-range order and the tran- 
sition temperature will depend on the energy scale as 
1/ In (eo/e), rising very steeply with increasing e.^^ Our 
Monte Carlo simulations show that unless the spin tran- 
sition is significantly suppressed by fluctuations the spin 
and orbital transitions would happen together. 

The simulation results also show that the isolated or- 
bital transition is in the universality class of the 2D Ising 
model. While we have not been able to observe the iso- 
lated finite temperature spin transition when it is sepa- 
rated from the orbital transition, on general grounds, we 
expect it also to be a continuous transition in the univer- 
sality class of the 2D Ising model (due to a small but non- 
zero spin space anisotropy). If the 3D couplings are more 
important than spin space anisotropy, then the transition 
could have a significant crossover region in the universal- 
ity class of the 3D Hcisenberg model (but will ultimately 
be in the 3D Ising universality class if a uniaxial spin 
anisotropy is also present). When the two transitions 
come together, the simulations find that the transitions 
tend to become first order. 

Based on the above, we propose a phenomenological 
phase diagram (See Fig. [10]) with the following features: 

1. The structural transition is driven by orbital order- 
ing, which happens at temperature Tq- It is set by the 
exchange energy scale in the problem. 



2. Let eo be the energy scale below which the long- 
wavelength fluctuations are suppressed. Then the ratio 
of Neel to orbital transition can be parametrized as (x — 



To 



2-x 
1 + In l/x 



and. 



Tn/Tc 



1 



for a; ^ 1, 



for a; 3> 1. 



(8) 



(9) 



We have two continuous phase transitions for x <^ 1 and 
one simultaneous first order transition for a; 3> 1. In 
between, the region a: ~ 1 can have a small stretch where 
the two transitions are practically inseparable but remain 
continuous (AB in Fig. [TU)) . The two transitions merge 
into a first order transition at the point B in Fig. (TU] 

We last note that in principle, doping can be the source 
of another kind of additional fiuctuation which signifi- 
cantly reduces both spin and orbital transition temper- 
atures from the mean-field values. This can also lead to 
the separation of structural and magnetic transition as 
observed in many families of iron pnictides. 



Discussion of materials 

We next discuss the relevance of this study to various 
experimental findings in the iron-pnictide materials. 

As mentioned previously, the parent compounds of the 
1111 family have two separate second order phase tran- 
sitions, while in the 122 family the two transitions are 
closer to each other in temperature. In the undoped 
BaFe2As2, the two transitions are slightly separated, 
where the structural transition starts as second order and 
is followed by a simultaneous first order jump both in the 
lattice distortion and magnetic transition at a lower tem- 
perature. On the other hand, in CaFe2As2 and 
SrFe2As2 the structural and magnetic phase transitions 
happen together as a single first order transition 
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The three behaviors reported in different iron-pnictide 
materials are all captured by our phase diagram of a cou- 
pled spin-orbital Hamiltonian. In particular, phase tran- 
sitions in the 1111 family correspond to the case when 
e is small and away from eo, where the structural and 
magnetic phase transitions are separated and of second 
order. On the other hand, phase transitions of CaFe2As2 
and SrFe2As2 correspond to the case when e is much 
larger than eo, where the two transitions occur as a sin- 
gle first order transition. The region e ~ eo is relevant 
to BaFe2As2. In this case, the structural and magnetic 
transition temperatures can be very close and there is a 
tendency for the magnetic transition to become first or- 
der. This indicates that these materials are close to the 
boundary between the distinct regions. 
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One can further ask which interaction term controls e 
in the iron-pnictide materials. In this study we have in- 
vestigated the role of spin space anisotropics described by 
Eq. © or Eq. ((T]). Their effects on the transition tem- 
perature are in essence the same as exchange couplings 
in the third direction. (38, 51, 53| Phenomenologically, e 
would refer to the larger of the terms in determining the 
phase transitions. It is known that the 122 family is more 
disperse in the third direction than the 1111 family. In 
particular, in the 122 family spin-wave spectra from neu- 
tron scattering are usually fitted with an additional 3D 
exchange coupling Jc, while for the 1111 materials Jc is 
essentially zero. In BaFe2As2, the third-direction cou- 
pling is non-zero but also appears small; the reported 
Jcl J\ is roughly 1%. 77, ll] Since a 2D Ising universal- 
ity has been found for this material, [t^ a uniaxial spin 
anisotropy could be more important. On the other hand, 
in CaFe2As2 and SrFe2As2, Jr is more substantial and 
Jcl J\ is estimated to be 10%. fi^,!?!] Therefore, in these 
materials spin exchange coupling in the third direction 
could be the controlling factor for e. 

We note that coupling to other degrees of freedom such 
as the lattice variable could also turn the isolated orbital 
or magnetic transition into first order. However, besides 
the phase diagram, the orbital variables have proven in- 
dispensable in describing various other properties of iron- 
based superconductors such as the emergent transport 
anisotropics. 
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A modest orbital polarization has been reported by 
angle- resolved photoemission (ARPES) experiments per- 
formed on the 122 family of iron pnictide materials. [331- 
IsgI ] This observation is crucial in explaining the strik- 
ing phenomenon that in these materials the resistivity 
is smaller in the longer AFM axis. This un- 

expected behavior is striking especially because optical 
measurements indicate a smaller scattering rate along 
the shortened FM direction, [i^, [l^ It is the presence 
of an anisotropic effective mass due to a preferred occu- 
pation of 

^xz over d,yz orbitals on the Fermi level that 
renders a better conducting pathway along the AFM 
direction. [H ill 

As mentioned previously, with a preferential occupa- 
tion of d:cz orbitals over dyz orbitals, relativistic spin- 
orbit coupling can induce an orbital angular momentum 
in the xy plane and lead to a single ion anisotropy. An 
excess population of d^z orbitals (through an induced 
dxz +idxy piece) can favor spins to point along the x-axis 
while excess population of dyz orbitals can favor spins to 
point along the y direction. We propose that this mecha- 
nism is the reason why the observed directions of ordered 
spin moments are tied to antiferromagnetism and end up 
point along the AFM direction. 

We last note that a possible orbital ordering has also 
been proposed for Fei+„Te:^Sei ^ (the so-called 11 fam- 
ily of iron chalcogenides) . [83|, [SJ]. In these materials, 
the ordered moments form a (7r/2,7r/2) diagonal double 



stripe pattern, and the spin orientation points toward the 
FM direction. 8^, 8^ Based on our discussion above, we 
believe this implies on the Fermi level a preferred pop- 
ulation of Wannier functions whose orbital lobes point 
along the same direction. One direct consequence of this 
prediction is that in 11 iron chalcogenides resistivity is 
smaller in the FM direction, [j^ . Ist} This is indeed consis- 
tent with recent resistivity measurements. The above 
prediction could be further tested by future ARPES and 
optical experiments on de-twinned iron chalcogenides. 



CONCLUSION 

In summary, we have studied finite-temperature phase 
transitions in a Hamiltonian of coupled Heisenberg spin 
and Ising orbital degrees of freedom. Using mean-field 
theory, Monte Carlo simulations, and general arguments 
we established the phase diagram of such a spin-orbital 
model and discussed its relevance to the iron-pnictide 
superconductors. We found that if spin rotational invari- 
ance is preserved, the magnetic transition temperature 
is pushed to zero in accord with the Mermin- Wagner 
theorem. In this case, there is only one single finite- 
temperature orbital phase transition which belongs to 
the 2D Ising universality class. By introducing spin space 
anisotropics into the Hamiltonian, spins can order at fi- 
nite temperatures and the magnetic and orbital transi- 
tions are found to couple together and become first or- 
der. This phase diagram captures several observed be- 
haviors in the 1111 and 122 families of iron pnictides. 
We also studied the case when relativistic spin-orbit cou- 
pling leads to a uniaxial anisotropy and found that the 
preferred spin-orientation is driven by orbital order. This 
explains why the direction of ordered moment in these 
materials is tied to their antiferromagnetism. 

In the field of iron-based superconductors, there are 
several open questions that remain to be answered. It 
is interesting to further explore other experimental im- 
plications of model Hamiltonians with coupled spin and 
orbital degrees of freedom. For example, can fluctua- 
tions in orbital and/or spin variables account for various 
anomalous phenomena that occur above the structural 
and magnetic transition temperatures? [89| What are the 
effects of orbital order and orbital fluctuations on twin 
boundaries? Are they related to the enhanced supercon- 



91| 



ductivity at domain walls in these materials? [& 
Calculations to address these interesting open questions 
are areas of future study. 
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FIG. 11: The acceptance probability is plotted versus temper- 
ature for several feedback steps. The initial geometric distri- 
bution shows a pronounced dip in acceptance near the critical 
point and later feedback steps show how this is corrected by 
clustering replicas around Tc- 



APPENDIX 

Monte Carlo Simulations with Parallel Tempering 

We have used a parallel tempering Exchange Monte 
Carlo (EMC) method to simulate our models. 
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It is an efficient extended ensemble simulation method 
that simulates multiple copies (replicas) of the system 
simultaneously at different temperatures. Exchanges be- 
tween replica configurations are accepted or rejected in 
accordance with detailed balance. Replica exchange has 
been used to study systems spanning many fields includ- 
ing strongly correlated systems, biological pathways, and 
spin glasses, [o^] The advantage of these methods is that 
while at high temperatures the system's memory is erased 
and when replicas go back to lower temperatures they 
explore large phase space uncorrelated in monte carlo 
time, [oil 

Recently, Katzgraber et al. [o^ showed that in order 
to maximally benefit from EMC, the temperature dis- 
tribution must be determined in a nontrivial way via a 
"feedback" method. The temperature distribution {Ti} 
is obtained by starting with some initial set and record- 
ing statistics on the "round trip" time from Tiow to Thigh- 
Minimization of this round trip time results in the opti- 
mal distribution. The endpoints, {Tiow -Thigh} are fixed 
and feedbacks of the simulation are done until the dis- 
tribution converges. The evolution of acceptance proba- 
bility is shown in Fig. 1111 Once {Ti} is determined, an 
exchange monte carlo simulation is performed using the 
stored optimal temperature distribution. 

Between the EMC moves that exchange replicas, one 



has some freedom in how to update each individual 
replica, provided one can always know the energy of that 
replica. We choose to do local spin/orbital flips by sweep- 
ing over the lattice and randomly choosing whether or not 
a site in the lattice attempts a spin flip or an orbital flip. 



Orbital and Spin Measurements 

There are two order parameters of interest for our 
Hamiltonian, one associated with the orbital degrees of 
freedom and another with the spin degrees of freedom. 
We measure second and fourth moments for both vari- 
ables, and, in the case of the magnetization, we measure 



at the two AFM wave vectors of interest, Qi 
and Q2 = (0,7r). 
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The fii take value or 1 in our model and (n) = ^ . The Si 
are classical Heisenberg spins with magnitude unity and 
{S) = 0. These orbital and spin measurements are used 
to evaluate Binder cumulant ratios as discussed below. 



Binder Ratios 

We deflne the Binder ratios for the orbitals <?„ and for 
the spins gs through the relations: 
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At low temperatures, the spin and orbital order param- 
eter distributions will be sharply peaked at their ex- 
tremum values. At high temperatures all variables will 
have gaussian distributions. For T ^ Tc the orbital 
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quantities are: 
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For the spins, the low temperature hmits are 
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The two orbital orders divide the system between Qi and 
Q2, resulting in different limits than a system without 
competing ordering wave vectors. 

At high temperature, T ^ T^, we get well known re- 
sults for the Binder ratio: 
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The difference between orbitals and spins comes purely 
from the dimensionality of the variable. We note that for 
gs, in contrast to the low temperature limits, the high 
temperature limits do not depend on the presence two 
ordering wavevectors. 



First and Second Order Phase Transitions 

We can estimate the thermodynamic Tc by carefully 
studying the size dependence of various physical quanti- 
ties. We rely on the Binder ratios defined previously and 
well known finite size scaling arguments to address the 
types of transitions we measure. We propose the usual 
scaling ansatz for the susceptibility. 



t = 



T-Tr_ 
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x{t.L)^Lix{L^\t\). 
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t is the reduced temperature, and x(^: L) is the suscepti- 
bihty per spin for a system of size L. % is some unknown 
but universal function and v and 7 are critical exponents 
which denote the power law divergence at Tc- 
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FIG. 12: For an anisotropic ( A = 0.75 ) 24 x 24 system, we 
plot the orbital and spin binder ratios, gn and gs- gn and 
gs both develop what seem like a divergence at the transition 
temperature. The development of such an incipient diver- 
gence is an indicator of a first order transition. 



The Binder ratios g„ and gs have the property that at 
Tc, they are independent of system size, universal con- 
stants of the system. We find Tc from Binder Ratio mea- 
surements for many pairs of systems of size L and 2L and 
plotting versus temperature. A crossing for a given pair 
gives a constant and Tc- Size dependence of this constant 
exponentially decays versus the system size. [9^ A more 
prominent size dependence occurs for the spins than or- 
bitals in Fig. [T] We extrapolate the size dependence to 
large L by fitting the exponential decay. The y-intercept 
of this fit is the thermodynamic Tc- 

Next we discuss the determination of the order of the 
transition that motivates our phase diagram for the pnic- 
tides. At a second order phase transition, various ther- 
modynamic quantities develop power law singularities 
characterized by critical exponents, in this case v and 
7. We arrive at Fig. 3] by varying critical exponents 
until a collapse of all points is achieved. In the case 
of anisotropy, there are no critical exponents that pro- 
duce a good data collapse, an indication that the transi- 
tion is not second order. To support the claim that the 
anisotropy leads to first order transitions, we show a plot 
of binder ratios for spin-space anisotropy (A = 0.75) in 
Fig |12l The binder ratio for spins develops a divergence 
near Tc that is accompanied by a weak divergence for the 
orbital binder ratio. On its own, this method does not 
conclusively establish the first order nature of the tran- 
sition. However, in conjunction with the lack of critical 
exponents, we propose the phase diagram in Fig. 1101 
Larger system sizes would be helpful in further stud ying 
the divergence of binder ratios in Fig|9]and Fig[T2l 97 [ 



